Computational Seismology
Reproducible Papers - Syngine Paper
In [1]:
%matplotlib inline
import obspy
import matplotlib.pyplot as plt
import numpy as np
plt.style.use("seaborn-whitegrid")
import copy
import io
import instaseis
import json
import requests
In [2]:
SYNGINE_URL = "http://service.iris.edu/irisws/syngine/1/query"
In [3]:
network = "IU"
station = "ANMO"
# Get station information from the IRIS FDSN service.
from obspy.clients.fdsn import Client
c = Client("IRIS")
print(c.get_stations(network=network, station=station, format="text")[0][0])
In [4]:
# The param file is only used to extract the source parameters. This is
# thus consistent with the other figures but can of course also be done
# differently.
filename = "chile_param.txt"
# Parse the finite source wiht instaseis.
finite_source = instaseis.FiniteSource.from_usgs_param_file(filename)
# Compute the centroid of it.
finite_source.compute_centroid()
# src is now the centroid of the finite source.
src = finite_source.CMT
# Common query parametersh su
params_common = {
# IU.ANMO
"receiverlatitude": 34.95,
"receiverlongitude": -106.46,
"dt": 0.1,
"origintime": src.origin_time,
"components": "Z",
"model": "ak135f_2s",
"format": "miniseed",
"units": "velocity"}
# Parameters only needed for the point source.
params_ps = copy.deepcopy(params_common)
params_ps["sourcelatitude"] = src.latitude
params_ps["sourcelongitude"] = src.longitude
params_ps["sourcedepthinmeters"] = src.depth_in_m
params_ps["sourcemomenttensor"] = ",".join(
str(getattr(src, _i)) for _i in ("m_rr", "m_tt", "m_pp", "m_rt", "m_rp", "m_tp"))
print(finite_source)
print(finite_source.CMT)
In [5]:
import copy
import collections
seis = collections.OrderedDict()
source_widths = [2.5, 5, 10, 25, 50, 100]
# Request one seismogram for each source with.
for sw in source_widths:
p = copy.deepcopy(params_ps)
# The sourcewidth parameter steers the width of the STF.
p["sourcewidth"] = sw
# Send it alongside.
r = requests.get(url=SYNGINE_URL, params=p)
assert r.ok, str(r.reason)
# Get the data and parse it as an ObsPy object.
with io.BytesIO(r.content) as f:
tr = obspy.read(f)[0]
seis[sw] = tr
# Plot only some phases.
tr.slice(tr.stats.starttime + 1000, tr.stats.starttime + 1500).plot()
In [6]:
import matplotlib.gridspec as gridspec
# Plotting setup.
fig = plt.figure(figsize=(10, 3))
gs1 = gridspec.GridSpec(1, 1, wspace=0, hspace=0, left=0.05,
right=0.62, bottom=0.14, top=0.99)
ax1 = fig.add_subplot(gs1[0])
gs2 = gridspec.GridSpec(1, 1, wspace=0, hspace=0, left=0.65,
right=0.94, bottom=0.14, top=0.99)
ax2 = fig.add_subplot(gs2[0])
plt.sca(ax1)
# Now plot all the seismograms.
for _i, (sw, tr) in enumerate(seis.items()):
tr.normalize()
plt.plot(tr.times(), 2.0 * tr.data - _i * 3, color="0.1")
plt.legend()
plt.xlim(0, 2000)
plt.yticks([0, -3, -6, -9, -12, -15], [str(_i) for _i in source_widths])
plt.ylim(-17, 2)
plt.xlabel("Time since event origin [sec]")
plt.ylabel("Source width [sec]")
plt.sca(ax2)
# Use an internal instaseis function to get the used STF.
from instaseis.server.util import get_gaussian_source_time_function
dt = 0.01
# Plot all the source time functions.
for _i, sw in enumerate(source_widths):
sr = get_gaussian_source_time_function(sw, dt)[1]
#sr = np.concatenate([sr2, np.zeros(1000)])
alpha = 0.4 - _i * 0.4 / len(source_widths)
plt.fill_between(np.arange(len(sr)) * dt - sw, sr, color="0.0", alpha=alpha, linewidth=0)
if sw == 25:
plt.plot(np.arange(len(sr)) * dt - sw, sr, color="0.0", lw=2)
ax2.annotate('25 sec', xy=(5, 0.07), xytext=(8, 0.10),
arrowprops=dict(facecolor='black', shrink=0.05))
plt.grid(True)
plt.xlim(-20, 20)
plt.ylim(-0.0005, 0.16)
plt.xticks([-10, 0, 10])
plt.yticks([0, 0.04, 0.08, 0.12])
plt.xlabel("Time [sec]")
plt.ylabel("Slip rate [m/sec]")
ax2.yaxis.tick_right()
ax2.yaxis.set_label_position("right")
ax2.yaxis.set_tick_params(length=2)
ax2.yaxis.set_tick_params(pad=4)
ax2.xaxis.set_tick_params(length=2)
ax2.xaxis.set_tick_params(pad=4)
ax2.xaxis.set_tick_params(color="#CCCCCC")
ax2.yaxis.set_tick_params(color="#CCCCCC")
plt.savefig("source_width.pdf")
